Nearsightedness of Electronic Matter in One Dimension 
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The concept of nearsightedeness of electronic matter (NEM) was introduced by W. Kohn in 1996 
as the physical principal underlining Yang's electronic structure alghoritm of divide and conquer. 
It describes the fact that, for fixed chemical potential, local electronic properties at a point r, like 
the density n(r), depend significantly on the external potential v only at nearby points. Beyond 
a distance R, changes Av of that potential, no matter how large, have limited effects on local 
electronic properties, which tend to zero as function of R. This remains true even if the changes in 
the external potential completely surrounds the point r. NEM can be quantitatively characterized 
by the nearsightedness range, R(r, An), denned as the smallest distance from r, beyond which any 
change of the external potential produces a density change, at r, smaller than a given An. The 
present paper gives a detailed analysis of NEM for periodic metals and insulators in ID and includes 
sharp, explicit estimates of the nearsightedness range. Since NEM involves arbitrary changes of the 
external potential, strong, even qualitative changes can occur in the system, such as the quantization 
of the energy bands or the filling of the insulating gap of an insulator with continuum spectrum. In 
spite of such drastic changes, we show that Av has only a limited effect on the density, which can 
be quantified in terms of simple parameters of the unperturbed system. 

PACS numbers: 71.10.-w, 71.15.-m 



I. INTRODUCTION 

This paper is based on a preliminary remark by W. 
Kohnji about a general concept called "nearsightedness 
of electronic matter (NEM)" and on a recent short report 
(PK), 2 which amplified that remark in various aspects 
but did not include detailed proofs. In the present paper, 
we select that part of PK dealing with non-interacting ID 
electrons and provide a full discussion, including detailed 
proofs. Future publications will amplify other sections of 
PK. 

By "electronic matter" we understand a system of 
many electrons with significant wavefunction overlap, in 
equilibrium under the action of a given external potential 
v(r). We shall consider the change of a local electronic 
property, like the electron density n(r), under the ac- 
tion of an arbitrarily strong potential perturbation w(r'), 
which vanishes inside a specified sphere, |r — r'\ = R. 
Note that we allow situations when the perturbation 
completely surrounds the point r. NEM states that the 
resulting density change at r, An(r, R), is bounded by a 
function An(r, R), 

An(r, R) < An(r, R), (1) 
independent of the amplitude or shape of w(r'), and that 



An(r, R) — > 0, monotonically as R — > oo, 



(2) 



provided only that the chemical potential /i is held fixed. 

The essence of NEM is contained in Eq. . Although 
it may not look very special, consider our perturbing po- 
tential w{r'), confined outside the sphere of radius R. 
The common sense says that, as we increase its strength, 
we need to increase the radius R if we want to maintain 
its effect at the center of the sphere below a certain level. 



The common sense also says that we need to increase R 
to infinity as we make w{r') stronger and stronger. In re- 
ality, this is not so: if the chemical potential is kept fixed, 
these effects will saturate and, in fact, no matter what 
w(r') we put outside the sphere of radius R, they cannot 
exceed a certain upper bound, which we will determine 
in this paper. 

For a given An, we can solve for R in An(r, R) = An 
and define the nearsightedness range R(r, An). The sig- 
nificance of R(r, An) is the following: any perturba- 
tion beyond R(r, An), of arbitrary shape and amplitude, 
cannot produce a density change at r larger than An. 
R(r, An) provides a simple, quantitative measure of near- 
sightedness. 

The above formulation of NEM often reminds people 
of Thomas-Fermi screening, sometimes even when we dis- 
cuss insulators. However, let us point a few facts. If 
one puts a charge inside the uniform charged electron 
gas and calculates the density response, he will find that 
the Thomas-Fermi exponential screening is valid only 
very near the impurity. Further away from the impu- 
rity, he will see Friedel oscillations, decaying as an in- 
verse power lawi^ These oscillations are not negligible. 
In fact, the Friedel oscillations were observed experimen- 
tally, thought not directly, one year after their theoretical 
prediction. 4 However, this was not realized until much 
later, when Walter Kohn made the connection between 
the two resultsA& He showed that there is a big discrep- 
ancy between the prediction of the Thomas- Fermi screen- 
ing theory and these experimental results and that a self- 
consistent calculation along Friedel's lines, of the density 
change due to impurities in coper, brings the theory and 
experiment to almost perfect agreement. The picture 
that emerged was that, in the asymptotic region, the 
screening only renormalizes the amplitude of the Friedel 
oscillations. The whole issue was considered very impor- 
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tant at that time, because it clearly demonstrated the 
existence of a sharp Fermi surface in real metals. We 
also like to mention one very well known fact in surface 
physics, where self-consistent calculations of metallic sur- 
faces showed that the effective potential goes to the bulk 
value extremely fast, typically within one or two layers. 
However, the density oscillations extend much further 
into the bulk and they can be viewed as the Friedel oscil- 
lations generated by the screened surface potential. With 
these being said, we hope that the reader will dissociate, 
right from the begenning, NEM from the Thomas-Fermi 
exponential screening and the nearsightedness range from 
the Thomas-Fermi screening length. 

Quantum gases display non-local density responses to 
local perturbations because of two factors. First, the 
effect of any local perturbation propagates to further 
distances through inter-particle interactions. This ef- 
fect can be regarded as classical, since it manifests, in 
the same way, in classical gasses. Secondly, there is a 
purely quantum effect, that steams directly from the un- 
certainity principle. This paper is concerned with this 
purely quantum effect, so it neglegts the inter-particle 
interation effects, entirely. The fact that NEM exists 
for non-interacting systems is extremely important. To 
understand why, let us go back in time and recall that, 
at the beginning of the electronic structure calculations, 
when the exact diagonalization was the method of choice, 
people were facing the so called "exponential wall" when 
trying to extend the calculations to larger systems: be- 
cause the number of operations in such calculations scales 
exponentially with the number of atoms, N , their appli- 
cability was, and is still limited to systems containg a 
few tens of atoms. Density Functional Theory (DFT)L£ 
provided a powerful alternative: because the number of 
operations in DFT calculations scales as N 3 , we can now 
solve the electronic structure for systems containing hun- 
dreds of atoms. However, electronic structure calcula- 
tions for biological and nano systems, or for extremely 
complex materials, involve thousands of atoms. At this 
scale, we start feeling the 'W 3 wall." Ab-initio quantum 
calculation for such complex systems will require a new 
generation of DFT algorithms, scaling linearly with the 
number of atoms. It is now generally accepted that NEM 
is the physical basis for these algorithms. 9 

The quest for linearly scaling algorithms was initiated 
by W.T. Yang, who was the first to argue that O(N) algo- 
rithms are possible^ The algorithm proposed by Yang is 
known by the name of Divide and Conquer (DC). There 
are now several reviews on the linear scaling electronic 
structure calculations. We will mention here the one by 
Goedecker— and the one by Wu and Jayanthi^ which, 
at the time of their publication, gave an exhaustive dis- 
cussion of O(N) methodologies. If we examine carefully 
these methodologies, they are all based on the same orig- 
inal idea, namely gluing together calculations done for 
smaller systems. What is different, is the representation 
and the way the size of these smaller systems (the trun- 
cation) is determined. For example, the real space ap- 



proaches will use the decay of the density matrix while 
the localized basis set approaches will use the overlap 
of these functions to judge how large these subsystems 
should be. 

Let us focus on the original implementation given by 
YangiiS Consider a self-consistent DFT iteration process, 
for a large quantum system. Each iterative step con- 
sists in calculating the density of a non-interacting elec- 
tron gas in equilibrium under a given effective potential 
(known from the previous iteration). In traditional ap- 
proaches, this requires a number of operations that scales 
as N 3 . DC algorithm, if it works as it is supposed to, re- 
quires a number of operations that scales linearly with 
N. It goes like this: The large system is devided in non- 
overlaping sub-regiones, which are then surrounded with 
buffer zones. A global chemical potential, /i, is fixed and 
the orbitals are calculated and populated up to /i, for 
each sub-region + buffer zone. The density in the buffer 
zones is discarded so, at this point, one has calculated a 
density for each sub-region and, by putting together all 
these sub-densities, one can construct the global density. 
The charge neutrality condition is then checked for this 
global density and the chemical potential is adjusted, if 
necessary. Note that charge neutrality must be satisfied 
by the entire system, not by each subsystem. In this way 
we have completed the DFT iteration step in a number 
of operations that scales linearly with the size of the sys- 
tem. Now, the question is how accurate is this algorithm? 
In fact, the most important question is, can we obtain 
arbitrary accuracy with this algorithm? To answer, we 
need to compare the density calculated for a sub-region + 
buffer zone and the density calculated for the entire sys- 
tem at once, at the same chemical potential /i. Now one 
can see why NEM can be regarded as the basis for this al- 
gorithm: the artificial termination at the outer boundary 
of the buffer zone, no matter how it is done, represents 
the change in the effective potential in our formulation of 
NEM. For example, such changes of the effective poten- 
tial occure when one calculates the density matrix and 
ignores the points outside a sub-region, or when one cal- 
culates the density and ignores the elements of a localized 
basis set that are centered outside a sub-region. Since we 
have no control on how the effective potential is modi- 
fied by such truncations, Eq. JJ| is paramount: it tells 
us that the effects of any artificial termination cannot 
exceed an upper bound. This upper bound is an intrin- 
sic characteristic of the system: is independent of the 
method of termination. Eq. tells us that, if we take 
the buffer zones large enough, the difference between the 
sub-density and the real density can be made smaller 
than any desired accuracy. 

Examples and estimates of the width of the buffer 
zones, together with a discussion of how to optimize this 
algorithm in 1, 2 and 3D and how the CPU time scales 
with the the desired accuracy can be found in Ref . We 
also like to mention that DC has been recently imple- 
mented to systems containing as many as 65,000 atoms 
and shown that it can generate electronic structures of 
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the same quality as a traditional approach will do for a 
system of, let us say, 10 atoms H The tests performed in 
this numerical work agree qualitatively with our theoret- 
ical predictions. 

There is another important issue related to DC. The 
ground energy in DFT is given by: 



E = Y^ej+E xc [n]- J v xc (r)n{r)-- J 



1 f n(r)n(r') 



(3) 



All the above terms can be calculated directly from the 
density, except the first one. However, this term is just 
the integral of the energy density, 



e(r)= ]T c#,-(r)| s 



(4) 



and we will show that e(r) is also nearsighted. As a conse- 
quence, within the DC algorithm, e(r) can be calculated 
with arbitrary precision, like the density. 

The goal of this paper is two-fold. We want to prove 
NEM, i.e. Eqs. Q and (J2J, for a simple system, which 
is the ID non-interacting, spin i fermions, in periodic 
potentials, and also want to show that one can ob- 
tain exact estimates of the nearsightedness range, which 
is extremely important for DC. The case of ID non- 
interacting electrons is important for several reasons. In- 
spite of its simplicity, it captures all the important as- 
pects of nearsightedness. This allows for a thorough in- 
vestigation of NEM, while keeping the technical aspects 
at a reasonable level. The ID non-interacting case is 
relevant for linear molecular chains when treated within 
DFT. 



II. THE STRATEGY 

We shall first develop general tools that will allow us 
to compute, for arbitrary perturbations, the asymptotic 
behavior of the density change, 



An(st) = 2 J2 hMs)l 



2 2£|tf(*)| a l 

e°<n 



(5) 



and Ge = {E — H) 1 denote the Green's functions of the 
unperturbed and perturbed systems, respectively, then 



An(x) 



1 



[G E -G° E ](x,x) dE, 



(6) 



where C is a contour in the complex energy plane, sur- 
rounding the occupied states. This can be seen from the 
eigenfunction expansions of G E and Ge and the residue 
theorem. Now, the eigenfunction expansions of G E and 
Ge are, again, not very useful. Instead, we will use the 
following compact representation: 



G E {x, x') = 



ip < (x < )ip > (x > ) 

w<,^>) : 



(7) 



where x< = min(x, x') and x> = max(a;,x'); ip<{x) and 
ip>(x) are the solutions of the Schrodinger equation at 
energy E, satisfying the boundary condition to the left 
and right, respectively, and W(ip<, ip>) is the Wronskian 
of the two solutions. For infinite systems, the case con- 
sidered in this paper, ip<(x) and ip>(x) are the solutions 
decaying at ±oo, respectively. We will always assume 
that E does not belong to the energy spectrum. When 
the contour C intersects the energy spectrum, such as for 
the case of metals, Ge will have a discontinuity at the 
point of intersection. Strictly speaking, this point must 
be excluded from C, which does not change the result 
of the integration. For all the other points of C, Ge is 
uniquely defined and given by Eq. Q. Later, we will use 
the reflection and transmission coefficients to construct 
extremely simple and compact expressions of the Green's 
functions (see Eqs. JTSJl, and JSTJl). 

The last step of our strategy will be to identify the 
special point in the complex energy plane that determines 
the asymptotic behavior of the integral Eq. ©. 

This strategy will require from us to go into the com- 
plex energy plane. We will make use of the analytic struc- 
ture of the Bloch functions and band energies derived in 
Ref. 0, which is briefly discussed in the next section. 
These analytic structure results can be generalized to 
linear molecular chains and even to 3D crystals?^ Also, 
the above expression for the Green's function, Eq. J7J|, 
can be generalized to linear molecular chains^ or to 3D 
crystals^ In fact, the entire strategy can be applied in 
2 and 3D, as it was already shown in Ref. 



where ip®(x), and ipi(x), 6i are the wave functions and 
the corresponding energies of the unperturbed and per- 
turbed systems, respectively. The factor 2 in front of the 
sums comes from the spin. The above expression is not 
very useful when dealing with the asymptotic behavior of 
An(x). Instead, we will work with an integral represen- 
tation. Why an integral representation? For answer, we 
point to the theory of special functions, where the func- 
tions are most often defined and introduced as infinite 
series but, with no exception, their asymptotic behavior 
is derived from equivalent integral representations. 

We can obtain an integral representation of An(x) by 
using the Green's functions. Indeed, if G E = (E — Hq)^ 1 



III. THE UNPERTURBED SYSTEM 

Throughout this paper, v(x) will be taken as a peri- 
odic, inversion-symmetric potential. Following is a brief 
discussion, largely taken from Ref. ITiL of the periodic 
Schrodinger equation. 

The solutions of the periodic Schrodinger equation, 

[-d 2 /dx 2 + v(x)]ip = Eip, v(x + b) = v(x), (8) 

are the well known Bloch functions ipk (k the wave vec- 
tor) which will be normalized as in Ref. Il4l Their funda- 
mental property is ipk (x + b) — e ipk (x) . When dealing 
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FIG. 1: The analytic domain of E n (X) and %p n> \(x) consists of 
the unit disk except a branch cut, represented by the dotted 
line. Panel a) refers to n even and panel b) to n odd. 



FIG. 2: Illustration of how E n (X) maps a) the complex A- 
plane into b) the complex E- plane (n even). 



with complex values of k, is much more convenient to 
work with the variable A = e lkb , instead of k. Thus, from 
now on, we will index the Bloch functions by A; their 
fundamental property becomes: 

V>a(z + &) = AV>a(z). (9) 

The parameter A relates to the energy of ip\ through the 
following equation: 

A 2 — 2fi(E)X + 1 = 0, (10) 

whith n(E) the Kramers' function^ By examining the 
fundamental property Eq. @ , one can see that the phys- 
ical states correspond to the case |A| = 1 (A on the unit 
circle), otherwise ip\(x) explodes either to x = oo or 
x = — oo. For A on the unit circle, the solutions behave 
like waves, thus it is appropriate to use the term Bloch 
waves. When discussing arbitrary values of A, however, 
it is more appropriate to use the term Bloch functions. 

The energy spectrum consists of energy bands, indexed 
here by n = 1, 2, . . ., which are separated by energy gaps. 
The energy bands can be computed by solving for E in 
Eq. HI U|l for all A on the unit circle. Due to the symme- 
try A — » 1/A in Eq. IjlOfl . we can and shall restrict A to 
|A| < 1, and view ij)\(x) and ipi/\(x) as two independent 
wavef unctions. For |A| < 1, it follows from the funda- 
mental property, Eq. ©, that ip\(x) decays to zero as 
x — > oo and ipi/\(x) decays to zero as x — ► — oo. 

If we restrict A to the unit disk, E uniquely determines 
A. The opposite is not true, instead E(X) is a multi- 
valued complex function, with branch points of order one 
at Ai, A2, .... Each A„ is real and < (— 1)™A„ < 1. 
The corresponding energy, E n = E(X n ), is also real and 
located in the n-th gap. Estimates of A n in the small gap 
and tight binding limits are given in Appendix A. E(X) 
can be represented on a Riemann surface with one sheet 
corresponding to one band. The n-th sheet can be taken 
as the entire unit disk, except a cut extending from A„_i 
to A n , as shown in Fig. 1. Near a branch point, E(X) 
behaves as the square root, 

E(X) = E n + 2a n (X/X n - I) 1 ' 2 + .... (11) 



The function E(X) on the n-th Riemann sheet will be 
denoted by E n {X). 

The integral in Eq. © will be mapped into the complex 
A-plane, by changing the variable from E to A. Thus, it 
will be important to understand how the contour C looks 
in this plane. It is more easy to understand how a given 
contour 7 in the A-plane is mapped in the complex energy 
plane, i.e. to construct the points E(X) when A sweeps 
through the points of 7. The mapping E n (X), from the 
unit disk to the complex £?-plane, is generic in one di- 
mension, in the sense that it is qualitatively independent 
of the periodic potential. The general picture is as fol- 
lows: E n (X) maps the unit disk into a domain T> n (see 
Fig. 2). The domains V n , n = 1, 2, . . ., are disjoint (with 
the exception of a possible common boundary) , and, all 
together, they cover the entire complex E'-plane. Now 
let us consider several contours. In Fig. 2, we used the 
same line-style for a contour and its image. The thick 
lines in Fig. 2b represent the energy bands. Now con- 
sider the contours C\ and C2 in the A-plane, starting and 
ending at zero and surrounding the branch cuts infinitely 
close. They are mapped into C[ and C' 2 , shown in Fig. 2b. 
This figure displays only a finite sector of the two con- 
tours, which extend from — ioo to +ioo. The unit circle is 
mapped in a loop that surrounds the n-th band, infinitely 
close. The segment from —1 to A„_i is mapped on the 
real axis, from the lower edge of the n-th band down to 
_E„_i. The segment from X n to +1 is also mapped on 
the real axis, from E n down to the upper edge of the n- 
th band. The domain T> n mentioned above, lies between 
the curves C[ and C' 2 . From the above information, one 
should be able to construct, qualitatively, the image of 
any other contour. 

ip\(x) and ipi/\(x) are multi- valued analytic functions 
of A, with branch points of order 3 at A„. For A on the 
n-th Riemann sheet, ipx( x ) wu l De denoted by ipn,\(x). 
Both functions diverge at the branch points as, 
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and 



(A/A„ - 1)1/4 



(13) 



where u n (x) and are periodic (antiperiodic) func- 

tions for n even (odd) and q„ is defined by |A„| = e~ q ™ h . 
The Wronskian of the two independent Bloch functions 
is given by 



bX dE n {\) 
2tt dX 



(14) 



Consequently, the Green's function G E satisfies the iden- 
tity: 

G E {x,x') — = 2i^ nl/x (x < )ip n> x(x > )^-. (15) 



IV. THE EFFECT OF PERTURBATIONS 
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FIG. 3: The contour of integration for metals in a) complex 
B-plane and b) complex A-plane. The contour in panel a) 
extends from — ioo to +ioo. 



and 



i/j n ,x(x) + R n (X)i(! nt i/x(x), x < -L 
T n (X)ipn,x(x), X > 0, 



(20) 



A. One Sided Perturbations 

We consider here perturbations that are either to the 
left or to the right of the point x, where we measure the 
density change An (a;). Let us assume that w is to the 
left of x. For convenience, we choose the origin of x at 
the right edge of w, so that w is confined in the interval 
[— L,0], with L arbitrarily large but finite. We calculate 
the particle and energy density changes at x > 0. 

As already mentioned, the density change An(x) is 
given by 

An(x) = — / [G E (x, x) - Gl(x, x)]dE : (16) 
m J c 

where G% and Ge are the unperturbed and perturbed 
Green's functions, respectively, and C is a contour in the 
complex energy plane, surrounding the eigenvalues below 
cf (see for example Fig. 3a). Similarly, the change of the 
energy density is 

Ae(x) = — f E[G E (x, x) - G° E {x, x))dE. (17) 
m J c 

We can focus on An(x) and give only the final results for 
Ae{x). 

We construct the perturbed Green's function from two 
independent solutions of the Schrodinger equation: 

[-d 2 /dx 2 + v(x) + w(x)]ip(x) = Eip(x). (18) 

Outside the interval \—L, 0], the solutions are linear com- 
binations of ipx an d ipi/x- As already mentioned, we need 
the solutions decaying to Too, which can be conveniently 
written in terms of the reflection and transmission coef- 
ficients: 



_ / T n (X)tp n s/ x (x), x < -L 



Vn ^ x{ ' 1 ^ n .i/x(x) + R+(X)^ n ,x{x), X > 



(19) 



where E was taken to be in T> n (Fig. 2). The Wronskian 
of the two independent solutions is 



^AT„(A) — . 



leading to the following useful identity, for x > 0: 

[G E - G E ](x,x)— = 2iR+(X)ip n ,x(x) 2 ^. 

m bX 



(21) 



(22) 



This identity, together with Eqs. (|16l) and l|17l) . shows 
that, for x > 0, An(x) and Ae(x) are completely deter- 
mined by the unperturbed wavefunctions and reflection 
coefficient, in a simple and universal way. 

Eq. (|22|l also provides the analytic structure of the 
reflection coefficient: -R^(A) has branch points of or- 
der 1 at A„_i and A„. If we go around these branch 
points, i?+ becomes R^-i and R„+i, respectively. In 
other words, R£ are different branches of a multi- valued 
function R + (X). Near the branch points, 



R+(X) = R+(X n ) + r+(A/A n - l) 1 ' 2 + 



(23) 



The poles of i? + (A), if any, are mapped by E(X) into the 
poles of Ge, i.e. the energies of the bound states. Since 
the bound states are located in the gaps, the poles of 
R + (X) are always located on the real axis and away from 
the branch cuts. 

By evaluating the Wronskian of ipx( x ) ana ^ ^i/x^ x ) ^ or 
x < — L and for x > and equating the two results, one 
can derive the following identity: 



T(A)T(1/A) + R + (X)R + (1/X) = 1. 



(24) 



For |A| = 1, this identity becomes |T(A)| 2 + |i?+(A)| 2 = 1, 
showing that |i? + (A)| < 1 for all A on the unit circle. 
Similar conclusion holds for R~(X). 
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FIG. 4: The exact (solid) and asymptotic (dashed) An(x) for 
the model of Eq. [vo = -2, Vo = 10, b = 1, first band 
20% filled]. 



1. Metals 

We assume the n-th band partially occupied and 
choose C in Eq. (|16|l as in Fig. 3, where X F is defined 
by e F = E n (X F ). We write 



x = y + mb, 



(25) 



with y restricted to the first unit cell. Mapping into the 
A-plane by using Eq. (|22J) and recalling the fundamental 
property of the Bloch functions, for x > 0, we obtain 



An(x) = j f R+{X)^ n , x {y) 2 X 
An integration by parts gives 
An(x) = 21m 



dX. 



(26) 



R+(X F )lpn,\ F (x) 2 



mb 
X 2m d 



b Jq, 2m dX 



[R£(X)iPnAy) 2 ]dX. (27) 



The integral is of order 1 jx 2 , as it can be seen from an- 
other integration by parts. We can conclude 



An{x) -» -Im[R+(X F )<ip nAF (x) 2 



for large x. Similarly, 



Ae(x) -» ^I m [i?+(A F )VvA F (x) 2 ]. 



(28) 



(29) 



Since |i?+(A_p)| < 1, the amplitudes of An(x) and Ae(x) 
cannot exceed, in the asymptotic limit, the upper bounds 



and 



An(x) -> -\ip n ,\ F (x)\' 



Ae(x) — > efAn(a;), 



(30) 
(31) 



FIG. 5: a) No bound states are present in the insulating gap. 
b) Bound states are present, but all the states below E n are 
occupied (solid lines) and all the states above E n are unoc- 
cupied (dashed lines) . c) States below E n are unoccupied, d) 
States above E„ are occupied. 

independent of the shape and amplitude of w(x). This 
proves NEM for metals and one sided perturbations. 

A comparison between Eq. I|28|l and an exact calcula- 
tion of An(x) for the perturbed Kronig-Penney modeU^ 



vt(x) 



Vo, x < 0, 



(32) 



is shown in Fig. 4. Note that the asymptotic regime 
starts from about two lattice constants away from the 
perturbation. 



2. Insulators 

We assume the first n bands completely filled. For in- 
sulators, the calculations arc more involved since the per- 
turbing potential w(x) may generate bound states in the 
insulating gap (= the gap above the n-th band). There 
can be a discrete or a continuum set of states inside the 
insulating gap. When the set is discrete, there are four 
distinct possibilities, as shown in Fig. 5. Let us analyze 
these cases first. 

No bound states in the insulating gap: We can take C[ 
(with opposite orientation, see Fig. 2b) as the contour 
of integration in Eq. (|16fl . Mapping into the complex A- 
plane and using the fundamental property of the Bloch 
functions, gives 



An{x) = X 



2/» 



2i 



RtW^nAyf 



2m- 



1 dX 

An 



(33) 



with C\ shown in Fig. 2a and y and m defined in Eq. 12511 . 
The integrand diverges at A„ as (A — An)" 1 / 2 but this sin- 
gularity is integrable. Away from the branch point, the 



FIG. 6: The exact (solid line) and asymptotic (dashed line) 
An(x) for the perturbed Kronig-Penney model Eq. 1321 [vo = 
—3, b = 1, Vo = 50, first band completely filled]. 



integrand is finite and (A/ A„) 2m becomes small as we fol- 
low the contour C\ towards A = 0. Thus, for large to, the 
main contribution to the integral comes from the region 
in the immediate vicinity of the branch point. Expanding 
the integrand near A n and keeping the leading term^S we 
find 



FIG. 7: Typical behavior of hx (solid lines) and h\/x dashed 
lines. E± denote the upper/lower edge of the insulating gap. 
K% can take values only in the shade areas. 



An(x) ^ R+{\ n )- 



dX 



Ci 



1 



u n (x) 



2 e -2q n x 



(34) 



dh 1/x /dE = -4> 1/x (0)- 2 A/x(x) 2 dx (38b) 



The integral is equal to -2B(2m, 1/2) (B 
tion) and behaves asymptotically as 
elude, 



Beta func- 
sJl-Kjm. We con- 



AnU)- -27?-(A„ )(g 



1/2 



u n (x) 



Similarly, 



Ae(x) -» -2R+{X n )E n 



2tt 



1/2 



Unix) e 



2„-2q n x 



(35) 



(36) 



i?+(A) = - 



h\/\ ipi/\(0) 



Mo) 



(37) 



As in Ref . ITil we choose the phase of the Bloch functions 
such that ip\(0) = "0i/a(O), so we can eliminate the last 
factor in the above expression. For E in the insulating 
gap, tp^(x), ipxix) and tpi/xi x ) are rea l functions and 
since 



dh</dE = -^<(oy 



ip^ (x) 2 dx 



(38a) 



dhx/dE = ^a(O)" 



ipx{x) 2 dx, 



(38c) 



it follows that h\ = —h\/x and and h\/x are decreas- 
ing functions of E. The typical behavior of hx and hi/\ 
is shown in Fig. 7. Now, if there are no bound states in 
the gap, hf and hx cannot be equal for any E in the gap. 
Then, since h\ increases while decreases with E, 
can take values only in the shaded area of Fig. 7, below 
0. Consequently, the right side of Eq. (|37|) is smaller or 
equal to 1, i.e. 



An implementation of Eq. (|35|l to the perturbed Kronig- 
Penney model described in Eq. I|32|) is given in Fig. 6. 
Note again that the asymptotic regime starts from one 
or two lattice sites from the perturbation. 

To end the proof of nearsightedness, we need to show 
that |i? + (A„)| cannot exceed an upper bound. Since the 
reflection coefficient is not evaluated on the unit circle, 
the inequality |i? + (A)| < 1 is no longer guaranteed. How- 
ever, if and hx denote the logarithmic derivatives at 
x = of ipf{x) defined in Eq. (|19fl . and oftp\(x), respec- and 
tively, then 



\R + W\ < i, 



(39) 



remains valid when E is in the insulating gap. 

We can then conclude that the amplitudes of An(x) 
and Ae(ir), in the asymptotic limit, cannot exceed the 
upper bound 



An(x) 



2n 



1/2 



=± ) u n (xfe- 2 ^ x 
xb 



Ae{x) -> E n An(x). 



(40) 



(41) 



This completes the proof of NEM for insulators and one 
sided perturbing potentials that do not generate bound 
states in the insulating gap. 

Bound states in the insulating gap: We show in Ap- 
pendix B that if 



w{x)tp±[x) 2 dx ^ 0, 



(42) 



where tp±{x) denotes the Bloch function at the up- 
per/lower edge of an energy band, then w generates 
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bound states above/below this band, even for infinitely 
small coupling constants. Thus, the presence of bound 
states in the gaps is not a rare occurrence in one dimen- 
sion. 

The asymptotic forms of An(x) and Ae(x) depend on 
how the bound states in the insulating gap are occu- 
pied. When all bound states below the branch point E n 
are occupied and the ones above E n are unoccupied, i.e. 
the situation illustrated in Fig. 5b, the asymptotic be- 
havior, Eqs. and lj5rj|) . of An(x) and Ae(x) remains 
unchanged. 

Consider now that there are unoccupied bound states 
below E n , as illustrated in Fig. 5c. Let ipo, of energy Eq, 
be such a state. For x > 0, ipo is equal, up to a fac- 
torization constant, to the exponentially decaying Bloch 
function of energy E (E = E n (X )): 



ifi (x) 



(1 



X 2 )A 



J \lpn,\ (x)\ 2 dx 



1/2 



where 



A = 



\Mx)\ 2 dx (A< 1). 



(43) 



(44) 



We have |Ao| = e~ q ° b , with qo strictly larger than zero, 
and ip\ {x) = e~ q ° x uo(x), with u (x + b) = (—l) n uo(x). 
Since go decreases as Eq moves away from E n , the first 
unoccupied state will have the slowest exponential decay, 
among all unoccupied states below E n . Thus, when the 
contribution of these states is subtracted from Eq. (|35|l . 
one finds that the asymptotic form of An(x) is deter- 
mined by the first unoccupied bound state: 



Ara(x) 



2(l-Ag)A 

Jo \i>n,X ( x )\ 2dx 



tto(x) e 



(45) 



with the index referring to the first unoccupied bound 
state in the insulating gap. Similarly, 



Ae(x) 



2(l-Ag)£ A 

Jo \^n,X ( X )\ 2dx 



u Q (x) 2 e- 2qoX . 



(46) 



Since A < 1, the amplitudes of An(x) and Ae(x) cannot 
exceed, in the asymptotic limit, the upper bounds 



An(x) 



2(1 -Ag) 

Jo \^n,\ ( x )\ 2 dx 



u Q (x) 2 e~ 2q ^, (47) 



and 



Ae(x) -> E An(x). 



(48) 



The results remain the same if, instead of unoccupied 
bound states below E n , there are occupied bound states 
above E ni as in Fig. 5d. In this case, the index will 
refer to the last occupied bound state. 

Continuum states in the insulating gap: We consider 
the case when w(x) fills the entire insulating gap with 



(a) 



n-th band 






L 




Re X 



1X1=1 



FIG. 8: a) The contour of integration in the complex _E-plane 
(Ci was introduced in Fig. 2). b) The same contour in the 
complex A-plane. 



continuum spectrum, such as when the insulator is in 
contact with an infinite metal. In this case, i?+(A) has 
a branch cut on the real axis. The states are considered 
occupied up to a Fermi energy, e^, which is in the in- 
sulating gap of the unperturbed insulator. We consider 
only the generic case when ep ^ E n and define Xf by 
€f = E„(Xf)- This Xf is located strictly inside the unit 
circle, as opposed to the case of metals. We also define 
q F > so that |A F | = e~ qFb . 

In Eq. (|16J) . we consider the contour of integration 
shown in Fig. 8a. Mapping into the complex A-plane 
and using again the fundamental property of the Bloch 
functions, 

An(x) = X 2 F m j J KmnAv? (J;)'™ 1 

(49) 

With the new variable q defined by X/Xp = e qb , we can 
write 

An(x) = -4A!™Im J R+(X + )^ x (y) 2 e- 2mqb dq, (50) 

g>0 

where A + = A + i0 + . The asymptotic behavior can be 
extracted from a simple integration by parts: 



An(x) -> -2Im[i?+(Aj 



_, . x1 y fl , A ,,(x) 2 



(51) 



By writing ipn,x F {x) as up{x)e qpX , with uf(x + b) 
(— l)™Wi?(x), we can conclude: 



An(x) -» -2lm[R+(X+)]u F (xf 
Similarly, 

Ae(x) -» -2lm[R+(X+)]e F u F (x) 2 '- 



-2qpx 



(52) 



(53) 



An implementation of Eq. (|52|l to the perturbed Kronig- 
Penney model Eq. I132H is shown in Fig. 9. Notice again 
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FIG. 9: The exact (solid) and asymptotic (dashed) An(x) for 
the model of Eq. [Vo = -5, v = -2, 6 = 1 and e F = 4]. 



that the asymptotic regime starts from two lattice con- 
stants away from the perturbation. 

To end the proof of NEM, we need to give an upper 
bound on the amplitudes of An(x) and Ae(x). The imag- 
inary part of -R^(Aj) is proportional to the local density 
of states, g(E : x), at E = cf and x — 0. Indeed, Eq. lt22"|) 
provides the following identity: 

Im[i?+(A+)]^„, AF ( 2 :) 2 = -L^Ini^+^x)], (54) 



2tt dqp 



leading to 



Im[fl+(A+)] = d " F ' d 1 f ' g{e F ,0). 



(55) 



Note that the coefficient in front of g(eF, 0) is determined 
by the unperturbed system. If we limit ourselves to the 
generic case of w's that generate finite densities of states 
at 6 f , then NEM follows from Eqs. (JUJ) and (55). For 
practical applications, we consider this argument suffi- 
cient. However, to achieve a full proof of NEM, we need 
to consider also the cases when g(E, 0) diverges (or be- 
comes extremely large) as E — > ep. For this special cases, 
the asymptotic form of Eq. (|5U[) cannot be extracted from 
a simple integration by parts and An{x) is no longer 
given by Eq. (|52[) : its specific functional form will de- 
pend on the type of singularity of g(E,0). This special 
situations will not be discussed here. 



B. Two Sided Perturbations 

We consider here the case when the point x, where we 
evaluate An{x) and Ae(x), has perturbing potentials wl 
to the left and wr to the right, as schematically shown 
in Fig. 10. For one sided perturbations, An{x) decays as 
x moves further and further away from the perturbation 
and, because of this simple picture, NEM is intuitive and 



simply to grasp. When left and right perturbing poten- 
tials are present, this simple picture is gone: there will 
be interference terms in An{x), whose amplitude remain 
constant in the region between the two perturbing poten- 
tials. In addition, wl + wr can induce a strong, qualita- 
tive change of the system, namely, the energy bands may 
become quantized. In spite all of these, we will show the 
following: for metals, the interference terms are not neg- 
ligible, but An(x) still remains bounded. For insulators, 
the interference terms are exponentially small and An(x) 
is given by the simple superposition of the left and right 
density changes. Similar conclusions hold for Ae(x). 

For convenience, we fix the origin in the middle of the 
interval that separates the two perturbing potentials and 
consider the distance, R, from the origin to the right/left 
edge of wl/r to be an integer of &, R = Nb. We are 
interested in the behavior, for x near the origin, of An(x) 
and Ae(x) when R — > oo. We follow our general strategy 
and derive first an expression for the Green's function on 
the interval [— R, R] . We look again for the solutions of 
the Schrodinger equation, 

[-d 2 /dx 2 + v(x) + w L (x) + wr(x)]i/j(x) = Eij)(x), (56) 

which decay at =Foo. On the intervals from — oo to the 
left edge of wl and from the right edge of wr to +oo, 
these solutions can be expressed as in Eqs. (11911 and Ij20|l . 
in terms of the total (corresponding to wl + wr) trans- 
mission and reflection coefficients. Then, one can use 
again the reflection and transmission coefficients to con- 
tinue these solutions inside the interval — R < x < R. On 
this interval, they take the following form: 



and 



~ Tr(X) 
T(A) 



[M*) + &rW1>i/\(?)] (57) 



W>i/x(x) + R£(\)il>x(x)], (58) 



where T LiR (X) and Rf R {\) are the transmission and re- 
flection coefficients of the left/right potentials, and T(X) 
is the total transmission coefficient, 



T(A) = 



f L (X)f R (X) 
l-R+(X)R+(\Y 



(59) 



All these coefficients depend on R: If R (X) denote the 
reflection coefficients when the right/left edge of w^/r 
is at x = (thus R% r(X) are independent of R), then 
Eq. J33 gives: 



R% R (X) = \™R% R {\) 



(60) 



The Wronskian of the two independent solutions is the 
same as given in Eq. (|21|l . From the two independent 
solutions, Eqs. (|57|l and l|58|l . and their Wronskian, we 
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R=Nb 


R=M> 




< ► 


< »> 

\ 








\ Wr / 



x = 



FIG. 10: The case when the point x (near the origin), where 
we evaluate An(x) and Ae(a:), has perturbing potentials to 
the left and to the right. 



and only if both \R^ R (X)\ are equal to 1. In this case, the 
energy band degenerates into discrete energy spectrum. 
Because the left hand side of Eq. (|rJ2)l is slowly varying 
compared to the right hand side, one can get the qualita- 
tive picture by setting the left hand side constant. If the 
amplitudes of R^ R (X) are equal to 1 on the whole unit 
circle, then Eq. (|62|l has 47V solutions, A^, distributed 
on the unit circle; the spacing between two consecutive 
solutions is 27r/47V + 0(1/7V 2 ). This is the picture in the 
A-plane. In the -E-plane, the discrete energies are given 
by E(\k). The spacing between two consecutive energies 
is 27r<9 A £(A fe )/47V + 0{\/N 2 ). 



derive, for — R < x < R, the following identity: 



G° E (x,x)}^ 



[G E (x,x) 

2iX 2N R+{X) 
1-X^R+(X)R R (X) 

2tX 2N R R {X) 
l-X™R+(X)R R (X) 
4iX iN Rj(X)R R (X) 
1-X^R+(X)R R (X) 



(61) 



<tpx(x)ipi/\(x)^—. 



We now can use Eqs. I|16|) and (|17|1 to find An(x) and 
Ae(x), for x near the origin. By using the fundamental 
property of the Bloch functions, we can understand the 
behavior of each term in the above identity. For A not on 
the unit circle, as it is the case in our integrals, the first 
term decay exponentially as x moves to the right; the 
second term decay exponentially as x moves to the left, 
but the third term is periodic, with period b. Fortunately, 
the amplitude of this term becomes smaller and smaller 
as the two perturbing potentials are moved apart. 

From Eq. (|filjl . one can easily obtain the energy spec- 
trum of the perturbed system. So, let us discuss first 
how the simultaneous presence of j R affects the energy 
spectrum. We are interested in the last occupied band 
(indexed by n), so, from now on, A is considered on the 
n-th Riemann sheet. The energies of the discrete states 
are given by the poles of the Green's function. From the 
identity Eq. (|61() , one can see that these poles correspond 
to those A satisfying the equation 



Rt(X)R R (X) = l/X 



IN 



(62) 



Since the discrete state energies are real, the solutions 
of Eq. (|62|l are always located on the unit circle or on 
the real axis, away from the branch cuts. Inside the unit 
disk we have |A| < 1; consequently, |l/A 4Ar | becomes very 
large in the limit R — > oo. Thus, if there are solutions of 
Eq. ijrj^ inside the unit disk, they must be located very 
close to the poles of either Rj^(X) or R R (X). In other 
words, these solutions are perturbations of the bound 
states generated by wl (or wr) alone, already discussed 
in the previous subsection. Poles on the unit circle exist if 



1. Metals 

The effects of band quantization will be the strongest 
for metals, since the Fermi energy lies inside the band. 
We calculate the density change by integrating Eq. I|61|) 
along the contour of integration shown in Fig. 3a and 
map the integral into the complex A-plane. The asymp- 
totic behavior of An(x), for large R, is determined by 
the behavior of the integrand near Xf and I/Xf (see 
Fig. 3b). In the immediate vicinity of these points, we 
can replace the slowly varying functions in the right hand 
side of Eq. 161f) [the reflection coefficients and the Bloch 
functions] with their value at Xf and 1/Xp- The integral 
then can be explicitly calculated and the result is: 



An(x, R) 



2 toh^fwg 

x [R+^ixf + R r ^ 1/Xf (x) 2 ] 

+ llm[Hl-X F N RiR- R )M XF {x)\ 2 , 



with the reflection coefficients evaluated at Xf- 

To prove NEM, we need to find an upper bound on 
the above expression. A simple analysis reveals that the 
largest density changes occur when \R L R \ = 1, i.e. when 
the band is quantized at the Fermi energy. In this case, 
we can rewrite Eq. I|63|) as: 



An(x, R) 



ilm 



tanlT 1 ^ {R\R R , 



l/2l 



(64) 



Rc[(R+R R y/ 2 ^ F ( 3 



+ -Im[ln(l-Ar4^)]|KWr. 

The density change, as function of R, has discontinuities 
every time when Xj^R^R^ = 1, i.e. when a discrete 
energy crosses the Fermi level. These discontinuities are 
finite: since |Im[tanh _1 (z)]| < tt/4 and |Im[ln(l - z)}\ < 
7r/2, for \z\ < 1, the amplitude of the asymptotic term of 
An(x, R) cannot exceed the upper bound 



An(x, R) 



2tt 



I^A F (x)| S 



(65) 
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independent of wt, m potentials and of the position of the 
Fermi energy relative to the discretized energies. Simi- 
larly, the amplitude of Ae(x, R) cannot exceed the upper 
bound 



Ae(x, R) ->■ e F An(x, R), 



(66) 



and this completes our discussion of NEM for metals. 

These upper bounds are optimal, in the sense that 
there are iol/r potentials (the worst scenario) that gen- 
erate a density and an energy density that are equal to 
these upper bounds. By comparing with the results of 
the previous Section, one can see that interference has 
non-trivial effects: these upper bonds are not simply the 
superposition of the left and right upper bounds. 



2. Insulators 

We consider first the situation when there are no bound 
states in the insulating gap. We can take C[ of Fig. 2 as 
the contour of integration in Eq (|16|l . which is mapped 
into Ci in the complex A-plane. For any A on this curve, 
the denominators in the right side of Eq. I|61|l goes ex- 
ponentially to 1 as R — + oo. Consequently, in this limit, 
the structure of the integrand becomes completely anal- 
ogous with the one studied in the previous subsection. 
The asymptotic behavior can be extracted as previously 
and the result is 



An(x, R) -> -2#+(A„ 



-2R R (\, 



2ir 

m 

2tt 

m 

1/2 



1/2 



u n (x) e 



2„-2 9ll R 



1/2 



u' n {x) e 



2„-2 9ll R 



-2R+(X n )R R (\ n ) (^) u n (x)u' n ( 



-4g„R 



(67) 



Thus, in the limit of large R, An(x, R) is just the sum 
of the independent density changes due to the left and 
right potentials, plus an exponentially small correction. 
From the previous subsection, we can conclude that the 
amplitude of An(x, R) cannot exceed, for large R, the 
upper bound 



An(x, R) -> 2 



1/2 



[xf 



[xf 



-2<}„R 



Similarly, 



Ae(x, R) E n An(x, R). 



(68) 



(69) 



For the case when there are bound or continuum states 
in the insulating gap, the conclusion is the same: for 
large R,, the density change near the origin is the sum 
of the independent changes induced by the left and right 
potentials. Upper bounds on An(x, R) can be trivially 
derived from the previous subsection. 



V. THE NEARSIGHTEDNESS RANGE 

The nearsightedness range R(x, An) was introduced as 
the range beyond which any perturbation, no matter how 
large, induces a density change at x less than the given 
An. The asymptotic R(x,An), in the limit An — > 0, 
can now be easily calculated from the upper bounds on 
An(x), derived in this paper. Since the periodic sys- 
tems are macroscopically homogeneous, the asymptotic 
R(x,An) will be independent of x. 

When solving for R in An(x, R) = Art, we first average 
An(x, R) over one unit cell. For metals, Eq. I|65|) leads to 
the following asymptotic expression: 



R(An) -> 1/An. 



(70) 



Such universal behavior is characteristic only to 1 dimen- 
sion; in higher dimensions, the nearsightedness range will 
depend on the average particle density^ 

For insulators and w's that generate no bound states 
in the insulating gap, Eq. (|68|l leads to 



R(An) -> — ln-^-, 
1 ' 2q„ An' 



(71) 



where 



4^2^" 



[u n {x) 2 + u' n {x) 2 } dx. (72) 



In the small gap and tight binding limits, h is completely 
determined by the exponential decay constant q n , n — > 
Aq n \/2/tt and n — > Ay/q n /irb, respectively. 

It is important to notice that the nearsightedness range 
does not depend on the details of the underlying potential 
v(x), but on some simple parameters that can be defined 
also for non periodic potentials. For example, q n can 
be identified with the exponential decay constant of the 
density matrix. 



VI. DISCUSSION 

The above analysis provides a quantitative analy- 
sis of the nearsightedness of electronic matter for non- 
interacting fermions, moving in one dimension under the 
action of periodic potentials. Although the simplest case 
possible, it allowed us to understand the different mech- 
anisms behind NEM. Although we cannot point to one 
simple and general physical explanation of NEM, it is 
now clear that NEM is due to a destructive interference 
not of the wave amplitudes but of density amplitudes rij 
associated with the single particle eigenstates i/jj. The 
asymptotic behavior of An(x) was found to be deter- 
mined by the reflection coefficient. For specific cases, 
the amplitude of An(x) cannot exceed an upper bound 
simply because, when evaluated at allowed energies, the 
reflection coefficients are always smaller than 1. More 
general, and now including 2 and 3D, one will find that 
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asymptotic behavior of An(x) is determined by certain 
elements of the scattering matrix and, for specific cases, 
the unitarity of the scattering matrix imposes certain up- 
per bounds. The situation is, however, more complicated 
when bound states appear in the insulating gap or when 
the bands become quantized. 

We have introduced a new concept, the nearsighted- 
ness range, R(x, An), which is well defined only because 
there is this upper bound on An(x). R(x, An) is a char- 
acteristic of the unperturbed system and gives a simple 
and effective measure of nearsightedness. For periodic 
metals, we found R(x, An) to have, in the asymptotic 
limit An — > 0, a universal expression, namely 1/ An. For 
insulators, R(x, An) is strongly dependent on the band 
structure but has a weak, logarithmic dependence on An. 

Although the estimates given in this paper can be ap- 
plied only to ID systems, we think we gain some knowl- 
edge that can be useful for more general situations. We 
are convinced that NEM exists in dimensions higher than 
1, where it can be quantified in a similar way. In par- 
ticular, we believe that a complete theoretical analysis 
and optimization of the O(N) divide and conquer algo- 
rithm is possible in all dimensions. Preliminary results 
in this direction have been already given in Ref. |2j. The 
one dimension analysis proved to be extremely useful by 
providing a viable strategy and some understanding of 
the effects of the bound states in the in insulating gap 
and of the band quantization on NEM. 



band. Since fj, n > 1, q n ~ \y/2(\ix n 
gether with Eq. (|A2j) . lead to 



1), which, to- 



y m * n (E- +1 -E+). 



(A3) 



We consider now a periodic potential V a (x — lb), 
where V a (x) vanishes for \x\ > c and has atomic levels 
E n = — n = 1, . . .. In the limit 6 — > oo, we show that 



1, SV^K 

— In 

b ebW n 



(A4) 



where W n is the width of the n-th energy band. For x in 
[—6/2,6/2] and \x\ > c, the solutions of the Schrodinger 
equation at an energy E = —k 2 are of the general form 

_ J a-(k)e~ kx + b-(k)e kx , x < -c 



^ i a+(k)e- kx +b+(k)e kx \ x > c, ^ 



with 



a+(k) \ _f (h , ( a-(fc) 
b + (k) J- J Wl 6_(fc) 



(A6) 



T(k) being the transfer matrix of the potential V a . The 
energy levels of V a correspond to the zeroes of T 22 (k), 
already denoted by k n . The Kramers function is given 
by 
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APPENDIX A 

We estimate here the exponential decay rate q n , re- 
lated to the branch point by |A„| = e~ 9nb . According to 
Ref.Q 



We estimate first the bandwidths. We look for the so- 
lutions of /i(fc) = ±1, which give the band edges. For 
6 large, the solutions of this equation must be located 
very close to the zeros of T 22 (k) since, otherwise, the sec- 
ond term in Eq. I|A7() becomes very large. We can then 
linearize, T 22 (k) ~ (k — k n )T2 2 (k n ) and neglect the ex- 
ponentially small term in Eq. (|A7|I . in which case the 
equation /i(fc) = ±1 can be trivially solved, leading to 



W n = 



8k n e 



-fc„b 



T 22 (k n )\ 



We now calculate E n — —k^, defined bjsi^ 



(A8) 



q n = - ln[|/x n | + y/fj,l - 1], 



(Al) 



where [i n is the Krammers function evaluated at the 
branch point E n , defined by d[i/dE\ E= g = 0. 

For small gaps, the behavior of n(E) inside the entire 
gap is well approximated by a quadratic function of E: 



tx(E) ~ (-1)" 



b 2 (E-E+)(E-E- +1 ) 
2(E~ +1 - E n ) 



(A2) 



where E^ is the upper/lower edge of the n-th band, and 
m* is the effective mass at the upper edge of the n-th 



Q^Ti 2 (k n ) ~-6T 22 (fc„). (A9) 



fc=fe„ 



For 6 large, the solutions of the above equation must also 
be close to the zeroes of T 22 (k). Linearizing T 22 (k), we 
find k n = k n — 1/6 and the Kramers function evaluated 
at kji is 



T 22 (k n )e kn 



2eb 



4fc„ 

ebW n 



(A10) 



Since [i n 3> 1, q n ~ i ln[2/i„] and Eq. I|A4(1 follows. 
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APPENDIX B 



where 



Let w(x) be a perturbing potential of finite support 
and such that 



2tt 



w(x)ip+(x) 2 dx > 



w(x)tp-(x) 2 dx < 0, 



(Bl) 



(B2) 



where tp±(x) is the Bloch function at the upper/lower 
edge of an energy band. We show here that, even for in- 
finitely small coupling constants, such potential will pull 
bound states out from the band. 

Let Hq denote the periodic Hamiltonian and H = Ho + 
71/7. It is known that H has a bound state at some energy 
E if and only if the operator— 



K E = ~/w 1/2 {E - Hq)' 1 ^ 



1/2 



(B3) 



has an eigenvalue equal to l£L Here, w 1 / 2 = w/lw] 1 ^ 2 . 
We show that, for any given energy E below/above the 
band, Ke has an eigenvalue equal to 1 for some positive 
7, which decreases to zero as E approaches the edges of 
the band, provided the condition Eq. I|B1|I / l|B2|) is satis- 
fied. 

If n is odd, the lower edge of the band corresponds 
to A = 1. We take an energy E below such band and 
let A, which is real and less than 1, be its corresponding 
A-parameter. Eq. I|15|) gives 

K E (x,x') = -^W^ ^^V ^I^ 2 , 

(B4) 

and we notice that dE/dX oc A — 1, for A — ► 1, i.e. the 
kernel Ke{x,x') diverges at A = 1. We can separate the 
diverging part by expanding 



^\-^{x<)^\{x- 



= ip-{x)il>-{x) 

+ {X-l)W x {x,x'). (B5) 



This provides the following decomposition, 

k E = 7 a(X)\i Pl )(ip 2 \+ 1 A(X), (B6) 



bXdE/dX'' 



(B7) 



^ S TS^« (B8) 



and 



(B9) 



fi(x) = w(x) 1 ' 2 ip-(x) 
f2(x) = | w(x)| 1 / 2- 0_(a:). 
The first term of Eq. (|B6jl diverges while the second one 
is analytic at A = 1. Now let 'J be given by 



# = ( 7 1(A) - 1)-Vi, (BIO) 
which is well defined for small 7. Then 

K E V = [l+ 7 a(A)(p 2 |( 7 i(A) - I)" (BH) 
In other words, Ke has an eigenvalue at +1, if 

l + 7a(A)(^ 2 |( 7 i(A) -irV) =0. (B12) 
We can rewrite this equation as 



(B13) 



If we denote the left side with F(j, A), then F(0, 1) = 
and <9 7 F(0, 1) — — J w(x)ip-(x) 2 dx > 0, i.e. the condi- 
tions of the analytic implicit function theorem are satis- 
fied, which means that, for any A near +1, there is always 
a solution 7(A) to the Eq. i|B12(l . Moreover, 



bXdE 
7(A) = ^dA 



w(x)ip- (x) 2 dx 



(B14) 



where the dots indicate terms of order o[(l — A) 2 ]. 7 is 
real and positive, for E below the band, and goes to zero 
as E approaches the band edge. 

The other possible cases, A = — 1 and n even, follow in 
the same way. 
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